Tarea 2: Frequentist Inference

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Fecha límite de entrega:

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Elaboración de Código

Objetivo

a la segunda tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la segunda parte del curso, los cuales se enfocan principalmente en inferencia estadística, diseño de experimentos y test de hipótesis. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Enlaces a videos de las clases:

Documentación:

Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 2 son las siguientes:

# Manipulación de estructuras
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# Para realizar plots
library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
# Manipulación de varios plots en una imagen.
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine

Elimine eval=FALSE para ejecutar las celdas.

Pregunta 1: Estimadores.

# Function to calculate bias
calculate_bias <- function(estimates, p_true = 0.5) {
  mean(estimates) - p_true
}

# sigma = 0
cte <- rep(0.5, 1000)
x <- 1:1000
data <- c()

for (n in x) {
  X <- rbinom(n, 1, 0.5)
  epsilon_sigma <- rnorm(1, 0, 0)
  p_hat <- epsilon_sigma + mean(X)
  data <- c(data, p_hat)
}
df_0 <- data.frame(x = x, p_hat = data)
bias_0 <- calculate_bias(data)

# Create individual plots
p1 <- ggplot(df_0, aes(x = x, y = p_hat)) +
  geom_line(color = "blue") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimator of p (sigma = 0)", x = "n", y = "Estimator of p") +
  theme_minimal()

# sigma = 1
cte <- rep(0.5,1000)
x <- 1:1000
data <- c()

for (n in x) {
  X <- rbinom(n, 1, 0.5)
  epsilon_sigma <- rnorm(1, 0, 1)  # sigma = 1
  p_hat <- epsilon_sigma + mean(X)
  data <- c(data, p_hat)
}
df_1 <- data.frame(x = x, p_hat = data)
bias_1 <- calculate_bias(data)

p2 <- ggplot(df_1, aes(x = x, y = p_hat)) +
  geom_line(color = "red") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimator of p (sigma = 1)", x = "n", y = "Estimator of p") +
  theme_minimal()

# sigma = 2
cte <- rep(0.5,1000)
x <- 1:1000
data <- c()

for (n in x) {
  X <- rbinom(n, 1, 0.5)
  epsilon_sigma <- rnorm(1, 0, 2)  # sigma = 2
  p_hat <- epsilon_sigma + mean(X)
  data <- c(data, p_hat)
}
df_2 <- data.frame(x = x, p_hat = data)
bias_2 <- calculate_bias(data)

p3 <- ggplot(df_2, aes(x = x, y = p_hat)) +
  geom_line(color = "green") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimator of p (sigma = 2)", x = "n", y = "Estimator of p") +
  theme_minimal()

# sigma = 4
cte <- rep(0.5,1000)
x <- 1:1000
data <- c()

for (n in x) {
  X <- rbinom(n, 1, 0.5)
  epsilon_sigma <- rnorm(1, 0, 4)  # sigma = 4
  p_hat <- epsilon_sigma + mean(X)
  data <- c(data, p_hat)
}
df_4 <- data.frame(x = x, p_hat = data)
bias_4 <- calculate_bias(data)

p4 <- ggplot(df_4, aes(x = x, y = p_hat)) +
  geom_line(color = "purple") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimator of p (sigma = 4)", x = "n", y = "Estimator of p") +
  theme_minimal()

# Print plots
print(p1)

print(p2)

print(p3)

print(p4)

# Print out the biases for each sigma
cat("Bias for sigma = 0:", bias_0, "\n")
## Bias for sigma = 0: 0.0003753956
cat("Bias for sigma = 1:", bias_1, "\n")
## Bias for sigma = 1: 0.0492002
cat("Bias for sigma = 2:", bias_2, "\n")
## Bias for sigma = 2: -0.0140122
cat("Bias for sigma = 4:", bias_4, "\n")
## Bias for sigma = 4: 0.01957655

Respuesta

Respuesta

Estos muestran alta variabilidad, tanto para altos como bajos valores de \(\sigma\). El único que mantiene cierta regularidad es \(\sigma = 0\). Esto se debe a que para \(\sigma = 0\) no hay ruido en el estimador, por lo que es esperable que este se acerque al valor real a medida que aumenta n. Vemos que los otros valores de sigma si se diferencian del valor real, aumentando esta diferencia promedio a medida que aumenta n (para \(\sigma = 1\) y \(\sigma = 2\) es cercana a \(\pm 4\), mientras que para \(\sigma = 4\) llega hasta \(\pm 15\)).

Pregunta 2: Intervalo de Confidencia

En las preguntas 2, 3 y 4 deberá trabajar con el dataset diabetes_prediction_dataset.

datos <- read.csv("diabetes_prediction_dataset.csv", header = TRUE)
head(datos)

Una forma de estimar la distribución de la media de una población es a través de la realización de la sampling distribution de una distribución cualquiera. El valor obtenido en cada sampling distribution nos entrega un estimador de la media que posee una determinada distribución de la población. Sabiendo esto, calcule la sampling distribution de las columnas bmi y blood_glucose_level, obteniendo el intervalo de confianza de \(95\%\) para cada una de las medias obtenidas desde la distribución. Para realizar este experimento considere los siguientes puntos:

Hints:

  • Para comparar comparar las 2 estimaciones puede ser útil normalizar bmi y blood_glucose_level a una misma escala.
  • Para realizar la sampling distribution podría serle útil el comando sample().
  • Puede ser útil generar un dataframe para graficar con ggplot2.

Respuesta:

# Definimos tamaños de muestreo
tamano_muestra = 100
n_muestras = 10000

# Inicializamos estructuras necesarias
list_mean_bmi = vector()
list_interval_bmi = vector()
list_typeCI_bmi = vector()
list_prop_bmi = vector()
sucess_bmi = 0

list_mean_bgl = vector()
list_interval_bgl = vector()
list_typeCI_bgl = vector()
list_prop_bgl = vector()
sucess_bgl = 0

# Obtenemos la media y desviación estándar de cada columna


# Sampling distribution, calculo del intervalo de confianza y proporción.
for (i in 1:n_muestras) {
  
}

# Generamos dataframes para plotear mas facilmente los datos.
# Plot de Intervalos de confianza
# Plot de Intervalos de confianza
# Plot de proporción de CI

Respuesta

Pregunta 3: Estimación de Máxima Verosimilitud

El objetivo de esta parte será obtener los parámetros de las distribuciones de la media y la mediana de bmi. Para responder la pregunta realice los siguientes puntos:

Cabe señalar que el método de máxima verosimilitud deberá ser programado por usted y no podrá utilizar librerías que entreguen el valor directo (por ejemplo la librería MASS).

Respuesta

# Obtenemos métricas de muestreo con repetición
means <- c()
medis <- c()
n_muestras <- 10000
tamano_muestra <- 100

for (i in 1:n_muestras) {
  
}

# dataframe con las medias y medianas de las muestras
# Histograma media
# Histograma mediana
# Media
# función log likelihood
ll_plot <- function(mu, sigma) {
  # ...
}

ll_plot <- Vectorize(ll_plot)

# definir espacio donde se va a evaluar ll_plot
mu_m <- 
sigma_m <- 
likelihood_m <- outer(X=mu_m, Y=sigma_m, ll_plot)

filled.contour(x=mu_m, y=sigma_m, z=likelihood_m, xlab=expression(mu), 
               ylab=expression(sigma), color = topo.colors)
# Mediana
# función log likelihood
ll_plot <- function(mu, sigma) {
  # ...
}

ll_plot <- Vectorize(ll_plot)

# definir espacio donde se va a evaluar ll_plot
mu_d <-
sigma_d <-
likelihood_d <- outer(X=mu_d, Y=sigma_d, ll_plot)

filled.contour(x=mu_d, y=sigma_d, z=likelihood_d, xlab=expression(mu), 
               ylab=expression(sigma), color = topo.colors)
likelihood_mean <- function(param) {
  # Definimos los parametros de entrada de la funcion
  mu <- param[1]
  sigma <- param[2]
  # Definimos la likelihood como la suma logaritmica de la función de densidad
  # ...
}

# Agrrgue el rango donde operará nlminb
nlminb(objective=likelihood_mean, start=c(, ), lower=c(, ), upper=c(, ))
likelihood_med <- function(param) {
  # Definimos los parametros de entrada de la funcion
  mu <- param[1]
  sigma <- param[2]
  # Definimos la likelihood como la suma logaritmica de la función de densidad
  # ...
}

# Agrrgue el rango donde operará nlminb
nlminb(objective=likelihood_med, start=c(, ), lower=c(, ), upper=c(, ))
# Muestra de medias


# Histogrmas
# Muestra de medianas


# Histogrmas

Pregunta 4: Bootstrap I

En esta pregunta se error e inetrvalo de confianza para la varianza de la columna HbA1c_level.

Las actividades por realizar son las siguientes:

Respuesta:

# Bootstrap
B <- 5000
largo <- length(datos$HbA1c_level)
output <- c()

for (it in 1:B) {
  
}
# Gráfico de puntos
# Histogrma
# Obtenemos error e intervalos de confianze
se_vars <-
z_a2 <-
alpha <-

# límite inferios
l.CI <- 
# límite superior
u.CI <- 

sprintf('El intervalo de confidencia de 95%% de las varianzas es: (%.3f,%.3f)', l.CI , u.CI)
sprintf('El SE de la varianzas es: (%.3f)', se_vars)

Respuesta

Pregunta 5: Bootstrap II

El siguiente código genera una regresión lineal de bmi con respecto a age, es decir, una función lineal de age, \(y(age) = b + m\cdot age\), que pretende determinar el valor de bmi.

linearMod <- lm(bmi ~ age, data=datos)
print(linearMod)
## 
## Call:
## lm(formula = bmi ~ age, data = datos)
## 
## Coefficients:
## (Intercept)          age  
##    23.15536      0.09945
m <- linearMod$coefficients["age"]
b <- linearMod$coefficients["(Intercept)"]

ggplot() +
  geom_point(aes(x=datos$age, y=datos$bmi)) +
  geom_line(aes(x=datos$age, y=datos$age*m+b), color="red") +
  ggtitle("Regresión lineal") +
  ylab("bmi") + 
  xlab("Age")

Repita el proceso de la pregunta 4 para estimar el error e intervalos de confianza de los parámetros de la regresión (\(m\) y \(b\)). ¿Qué indican los resultados de Bootstrap?¿Un valor bajo en el error significa que la regresión es buena?

# Bootstrap
B <- 500
largo <- 1000
output_m <- c()
output_b <- c()

for (it in 1:B) {
  
}
# Obtenemos error e intervalos de confianze
se_m <- 
z_a2 <-
alpha <-

# límite inferios
l.CI <- 
# límite superior
u.CI <- 

sprintf('El intervalo de confidencia de 95%% de las varianzas es: (%.3f,%.3f)', l.CI , u.CI)
sprintf('El SE de la varianzas es: (%.3f)', se_m)
# Obtenemos error e intervalos de confianze
se_b <- 
z_a2 <-
alpha <-

# límite inferios
l.CI <- 
# límite superior
u.CI <- 

sprintf('El intervalo de confidencia de 95%% de los sesgos es: (%.3f,%.3f)', l.CI , u.CI)
sprintf('El SE del sesgo es: (%.3f)', se_b)

Respuesta

 

A work by CC6104